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Abstract 

The concept of topological sensitivity (TS) is extended to enable simultaneous 3D reconstruction of 
fractures with unknown boundary condition and characterization of their interface by way of elastic waves. 
Interactions between the two surfaces of a fracture, due to e.g. presence of asperities, fluid, or proppant, are 
described via the Schoenberg’s linear slip model. The proposed TS sensing platform is formulated in the 
frequency domain, and entails point-wise interrogation of the subsurface volume by infinitesimal fissures 
endowed with interfacial stiffness. For completeness, the featured elastic polarization tensor - central to 
the TS formula - is mathematically described in terms of the shear and normal specific stiffness (/^s, Kn) 
of a vanishing fracture. Simulations demonstrate that, irrespective of the contact condition between 
the faces of a hidden fracture, the TS (used as a waveform imaging tool) is capable of reconstructing its 
geometry and identifying the normal vector to the fracture surface without iterations. On the basis of such 
geometrical information, it is further shown via asymptotic analysis - assuming “low frequency” elastic- 
wave illumination, that by certain choices of (/^s,/^^) characterizing the trial (infinitesimal) fracture, the 
ratio between the shear and normal specific stiffness along the surface of a nearly-panar (finite) fracture 
can be qualitatively identified. This, in turn, provides a valuable insight into the interfacial condition 
of a fracture at virtually no surcharge - beyond the computational effort required for its imaging. The 
proposed developments are integrated into a computational platform based on a regularized boundary 
integral equation (BIE) method for 3D elastodynamics, and illustrated via a set of numerical experiments. 
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1. Introduction 


To date, inverse obstacle scattering remains a vibrant subject of interdisciplinary research with applica¬ 
tions to many areas of science and engineering [36]. Its purpose is to recover the geometric as well as 
physical properties of unknown heterogeneities embedded in a medium from the remote observations of 
thereby scattered waveforms. Such goal is pursued by studying the nonlinear and possibly non-unique 
relationship between the scattered field produced by a hidden object, e.g. fracture, and its characteristics. 

From the mathematical viewpoint, the fracture reconstruction problem was initiated in [31] where, 
from the knowledge of the far-held scattered waveforms, the shape of an open arc was identihed via the 
Newton’s method. This work was followed by a suite of non-iterative reconstruction approaches such as 
the factorization method [e.g. ng, the linear sampling method (LSM) [29] and the concept of topological 
sensitivity (TS) [25] that are capable of retrieving the shape, location, and the size of buried fractures. 
Recently, a TS-related approach has also been proposed for the reconstruction of a collection of small 
cracks in elasticity [5]. 

A non-iterative approach to inverse scattering which motivates the present study is that of TS Esun]. 
In short, the TS quantifies the leading-order perturbation of a given misfit functional due to the nucleation 
of an infinitesimal scatterer at a sampling point in the reference (say intact) domain. The resulting TS 
distribution is then used as anomaly indieator by equating the support of its most pronounced negative 
values with that of a hidden scatterer. The strength of the method lies in providing a computationally 
efficient way of reconstructing distinct inner heterogeneities without the need for prior information on 
their geometry. Recently, m demonstrated the ability of TS to image traetion-free eraeks. Motivated by 
the reported capability of TS to not only image - but also characterize - elastic inclusions this study 
aims to explore the potential of TS for simultaneous imaging and interfacial characterization of fractures 
with contact condition due to e.g. the presence of asperities, fluid, or proppant at their interface. 

Existing studies on the sensing of obstacles with unknown contact condition reflect two principal 
concerns, namely: i) the effect of such lack of information on the quality of geometric reconstruction, and 
ii) the retrieval - preferably in a non-iterative way - of the key physical characteristics of such contact. 
The former aspect is of paramount importance in imaging stress corrosion fractures [28], where the 
crack extent may be underestimated due to interactions at its interface, leading to a catastrophic failure. 
To help address such problem, m developed the factorization method for the shape reconstruction of 
acoustic impedanee cracks. Studies deploying the LSM as the reconstruction tool [e.g. [161 [18], fhe other 
hand, show that the LSM is successful in imaging obstacles and fractures regardless of their boundary 
condition. As to the second concern, a variational method was proposed in [20] to determine the essential 
supremum of electrical impedance at the boundary of partially-coated obstacles. By building on this 
approach, devised an iterative algorithm for the identification of surface properties of obstacles from 
acoustic and electromagnetic data. Recently, [32] proposed a Fourier-based algorithm using reverse-time 
migration and wavefleld extrapolation to retrieve the location, dip and heterogeneous compliance of an 
elastic interface under the premise of a) one-way seismic wavefleld, and b) absence of evanescent waves 
along the interface. 

Considering the small-amplitude elastic waves that are typically used for seismic imaging and non¬ 
destructive material evaluation, the Schoenberg’s linear slip model [41] is widely considered as an adequate 
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tool to describe the contact condition between the faces of a fracture. This framework can be interpreted 
as a linearization of the interfacial behavior about the elastostatic equilibrium state [38] prior to elastic- 
wave excitation, which gives rise to linear (normal and shear) specific stiffnesses and kg. Here it is 
worth noting that strong correlations are reported in the literature munii] between {kg^kn) and surface 
roughness, residual stress, fluid viscosity (if present at the interface), intact material properties, fracture 
connectivity, and excitation frequency. In this vein, remote sensing of the specific stiffness ratio kgjk^ has 
recently come under the spotlight in hydraulic fracturing, petroleum migration, and Earth’s Critical Zone 
studies isniE]. By way of laboratory experiments [HEZIIT], it is specifically shown that kgjkn - often 
approximated as either one (dry contact) or zero (isolated fluid-filled fracture) - can deviate significantly 
from such canonical estimates, having fundamental ramifications on the analysis of the effective moduli 
and wave propagation in fractured media. A recent study [HHH] on the production from the Cotton Valley 
tight gas reservoir, using shear-wave splitting data, further highlights the importance of monitoring ksjk^ 
during hydraulic fracturing via the observations that: i) the correlation between proppant introduction 
and dramatic increase in kgjk^ can be used as a tool to directly image the proppant injection process; 
ii) the ratio kgjkn provides a means to discriminate between newly created, old mineralized and proppant- 
filled fractures, and iii) kg/kn may be used to monitor the evolving hydraulic conductivity of an induced 
fracture network and subsequently assess the success of drilling and stimulation strategies. 

In what follows, the TS sensing platform is developed for the inverse scattering of time-harmonic elastic 
waves by fractures with unknown geometry and contact condition in On postulating the nucleation 
of an infinitesimal penny-shaped fracture with constant (normal and shear) interfacial stiffnesses at a 
sampling point, the TS formula and affiliated elastic polarization tensor are calculated and expressed in 
closed form. Simulations demonstrate that, irrespective of the contact condition between the faces of 
a hidden fracture, the TS is capable of reconstructing its geometry and identifying the normal vector 
to the fracture surface without iterations. Assuming illumination by long wavelengths, it is further 
shown that the TS is capable (with only a minimal amount of additional computation) of qualitatively 
characterizing the ratio kg/k^ along the surface of nearly-planar fractures. The proposed developments 
are integrated into a computational platform based on a regularized boundary integral equation (BIE) 
method for 3D elastodynamics. Eor completeness, the simulations also include preliminary results on the 
“high”-frequency TS sensing of fractures with specific stiffness, which may motivate further studies in 
this direction. 

2. Preliminaries 

Consider the scattering of time-harmonic elastic waves by a smooth fracture surface T C C (see 
Eig. with a linear, but otherwise generic, contact condition between its faces T^. Eor instance the 
fracture may be partially closed (due to surface asperities), fluid-filled, or traction free. Here, is a ball 
of radius Ri - containing the sampling region i.e. the search domain for hidden fractures. The action of 
an incident plane wave on T results in the scattered field u - observed in the form of the total field 

u(.i) = u%) + m( 4), 4 e (1) 
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Figure 1: Illumination of a hidden fracture F G with specific stiffness K by plane waves. Thus induced wavefield is 
monitored over . 


over a closed measurement surface 5'°^® = 9^2, where ^2 is a ball of radius R 2 ^ Ri centered at the 
origin. The reference i.e. “background” medium is assumed to be elastic, homogeneous, and isotropic 
with mass density p, shear modulus p, and Poisson’s ratio v. 

Dimensional platform. For simplicity, all quantities in the sequel are rendered dimensionless by tak¬ 
ing p, p, and Ri as the characteristic mass density, elastic modulus, and length, respectively. 

Sensory data. In what follows, the time-dependent factor e^^^ will be made implicit, where uj denotes the 
frequency of excitation. With such premise, the incident wavefield can be written as u'^{^) = 
where k = uj/c signifies the wavenumber; c is the relevant (compressional or shear) wave speed; b E Cl 
is the polarization vector, and d e Cl specifies the direction of propagation of the incident plane wave, 
noting that Cl stands for a unit sphere. For each incident plane wave specified via pair {b^d), values of 
the total field u{^) are collected over 

Governing equations. With the above assumptions in place, the scattered field u{^) can be shown to 
satisfy the field equation and interfacial condition 


t±(C) = T H(4) - C e r±, 


( 2 ) 


complemented by the Kupradze radiation conditions [T] at infinity. Here In] = |n| = n+ — u~ signifies 
the crack opening displacement (COD) on F; = n^-C: where is the unit normal on F^ 

(see Fig. Q; K{^) is a symmetric, positive-definite matrix of the specific stiffness coefficients; tf^ = 
n^- C'.Vu^ denotes the free-field traction on F^, and C is the fourth-order elasticity tensor 


C = 2p 


jsym 


l-2n' 
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in which I 2 and stand respectively for the second-order and symmetric fourth-order identity tensors. 
Following the usual convention [13], the unsigned tractions and normals on a generic surface S (e.g. tf^n) 
are referred to S~ and affiliated normal n~ where applicable. 

Here it is noted that X, which accounts for the interaction between r+ and F” due to e.g. surface 
asperities, fluid, or proppant at the fracture interface, may exhibit arbitrary spatial variations along F in 
terms of its normal and shear components. In light of the fact that the primary focus of this work is “low” 
frequency sensing where the illuminating wavelength exceeds most (if not all) characteristic length scales 
of a fracture - de facto resulting in the spatial averaging of its properties, it is for simplicity assumed 
that i) the normal specific stiffness is constant along F, and ii) the shear specific stiffness is both constant 
and isotropic @111101139]. More specifically, it is hereon assumed that 

= ks{eg<S)eg) + kn{n0n), Id = 1,2, (3) 

where kg = const, and kn = const, are the respective specific stiffnesses in the tangential (eg(^)) and 
normal (n(^)) directions of F; \eg\ = 1; (8) signifies the tensor product, and Einstein summaton convention 
is assumed over repeated indexes. 

Cost functional. For the purposes of solving the inverse problem the cost functional is, assuming given 
incident wavefield defined as 

J(r.HaO = [ (4) 

J Sohs 

in terms of the least-squares misfit density 

C) = i . W{^) ■ , (5) 

where are the observations of nj^obs (say polluted by noise); v is the simulation of u computed for 
trial fracture FtHai: and FF is a suitable (positive definite) weighting matrix, e.g. data covariance operator. 

3. Topological sensitivity for a fracture with specific stiffness 

On recalling Q and denoting 

F, = c + ^r,,ai, (6) 

where Ftnai contains the origin, the topological sensitivity (TS) of the featured cost functional can be 
defined as the leading-order term in the expansion of J{Ts) with respect to the vanishing (trial) fracture 
size, 5^0 uni. In what follows, Ftnai is taken as a penny-shaped fracture of unit radius with unit 
normal n', shear specific stiffness and normal specific stiffness Kn- To maintain the self-similarity 
of Fg - in terms of mechanical behavior - with respect to its vanishing size, is further endowed with 
an 5-dependent matrix of specific stiffness coefficients, namely 

Ke = + ^(n'0n'), (7) 
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where (e']^,e 2 ,n') constitute an orthonormal basis. Considering the mechanics of rough surfaces in con¬ 
tact, one may interpret 0 in physical terms by noting that the height of asperities implicit to FtHai is 
downscaled by the factor of 5 when computing via 0 - see, e.g., Ch. 3 of m for a phenomenological 
derivation of the specific contact stiffnesses from the characteristics of a rough interface. 

On the basis of the above considerations, the topological sensitivity T(^°; n'^Kn^i^s) is obtained from 
the expansion 

J(r,) = J(0) +/(e)T(r;n',K„,K,) + o(/(£)) as e ^ 0, (8) 

where f{e) 0 with diminishing 5 , see also [4^ |23l US] [M] . Thanks to the fact that the trial scattered 
field v{^) = v{^) — u^(^) due to F^ vanishes as 5 ^ 0, Q can be conveniently expanded in terms of v 
about u^, see m- As a result Q can be rewritten, to the leading order, as 

J(r,)-J(0) ~ = /(£)T(r;n',K„,«,). (9) 

J^obs OV 

Adjoint field approach. At this point, one may either differentiate ([^ at t; = n* and seek the asymptotic 
behavior of i;(^) over or follow the adjoint field approach [e.g.HSlIH] which transforms the domain 
of integration in ^ from 5'°^® to F^ - and leads to a compact representation of the TS formula. The 
essence of the latter method, adopted in this study, is to interpret the integral in 0 through Grafh’s 
reciprocity identity [1] between the trial scattered field and the so-called adjoint field 'u(^), whose 
governing equations read 


V 


V-[C:Vvm+pajH{^)=0, ^ e 


u: 


V-[C:Vu]{O+puJ^u{^)=0, Ce 




^ G 5'°’^" 


( 10 ) 

subject to the Kupradze radiation condition at infinity. Here t and i denote respectively the adjoint- and 
scattered-field tractions; [u] = v~^ — v~ is the crack opening displacement on Fg, and 


[<1(0= linin(0-C'-(VM(^-7?n)-V-u(^-Fr?n)), $ e 8°^'' 


denotes the jump in adjoint-field tractions across 5'°’^® with outward normal n. Note that the adjoint field 
is defined over the intact reference domain, whereby u is continuous G and consequently 
on F^=^. As a result, application of the reciprocity identity over M^\F£ can be shown to reduce ^ to 

T(r;n',= {f{e)y^ [ (11) 

dr. 


3.1. Asymptotic analysis 

Considering the trial scattered field 'h(C), the leading-order contribution of the COD is sought on the 
boundary of the vanishing crack G F^) as 5 ^ 0. For problems involving kinematic discontinuities such 
as that investigated here, it is convenient to deploy the traction BIE framework m as the basis for the 
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asymptotic analysis, namely 


= n' C -.-f 


cl5^ 


puP'n'-C: [ J7(^,a;,a;)-(|i;]](8)n')(a:)dS'a;, 
-'r. 


( 12 ) 




where 


U = Ui{^,x,Lo)ek 0 Cj 


S = sf • (^, ®, w) efe 0 e* 0 ej ; 


U^{^,x,Lo) and T,^j{^,x,x) (given in Appendix A) denote respectively the elastodynamic displacement 


and stress fundamental solution due to point force acting at cc € in direction fc; f signifies the 
Cauchy-principal-value integral, and D^. is the tangential differential operator m on Fg, given by 

Dx{f) = -D/cz(/m) ^ ^ ^ki -^/cz(/m) ~ 

such that n'j^ = n'^{x) and fm,k = dfm/dxk- 


Scaling considerations. To facilitate the analysis of (12), all quantities are hereon assumed to be dimen¬ 


sionless. This is accomplished by taking the radius Ri of ball containing the sampling region (see 
Fig. [^, the mass density p, and the shear modulus p as the reference scales for length, mass density, and 
stress. Next, a change of variable cc = is introduced on T^ as motivated by ([^, which results in 

the scaling relations 


d5, = 5^d5^, tf{^) = tf{e)^ 0 {e), | G T,, 


(13) 


when 5^0. In this setting, the elastodynamic fundamental tensors in (12) are known to have the 
asymptotic behavior 


5](|,x,u;) = 5-^5 ](|,x) + 0(1), U{^,x,uj) = t/(|, x) + 0(1), x,^eT 


(14) 


in which U{^,x) and I](^,cc) signify the displacement and stress tensors associated with the Kelvin’s 
elastostatic fundamental solution m- Following the logic of earlier works [e.g.HOllillSl, the asymptotic 
behavior of a vanishing scattered field 


N(0 := I^Kr + eC) 


(15) 


as 5 ^ 0 can be exposed by substituting (13)-(14) into (12) which yields 

= e-^n'cU i:[l,x) : D^m{x) dS^ + 0 ( 5 ), | G F.^ai, (16) 

^rtrial 
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and seeking the balance of the featured leading terms [33]. To solve (16), consider a representation of the 
fracture opening displacement as 


[t^K^) C. e-4.(r)ivp(e, (17) 

where a{j are the components of the free-held stress tensor aj = C :Vu^, and {i,j = 1,2,3) are 


canonical solutions to be determined. On recalling that the free-field traction tf = n' ■ af in (16) is 


independent of 5, one immediately finds that a = 1 in (|17[) which reduces (16) to 


+ =n' C :j *):£>* [VP(S)d^*, | G T, 


(18) 


By analogy to the BIE formulation for an exterior (traction-free) crack problem in elastostatic [e.g. [TTl [8] 
one recognizes that, for given pair (i, j), integral equation ( jls] ) governs the fracture opening displacement 
due to tractions n'’{ci 0 ej -h ej O e^) applied to the faces of a “unit” fracture TtHai with 
interfacial stiffness iTtnai in an infinite elastic solid (recall that n'^ = Owing to the symmetry of 


with respect to i and j, (18) can accordingly be affiliated with six canonical elastostatic problems 


m 


The TS formula. Having {vj characterized to the leading order, one finds from @ @ and @ 
that /(e) = e^ and consequently 

T{C;n',Kn,Ks) = (Tf{^°):A:&{^°), A = ei®ej®( j |Vf-^(S)dS'*) ®n', (19) 

V ^Ftrial / 

where a = C \ Vii denotes the adjoint-field stress tensor, and A is the so-called polarization tensor - 
independent of and uj - whose evaluation is examined next. 


3.2. Elastie polarization tensor 

In prior works on the topological sensitivity [e.g.jSlIMliniES], relevant polarization tensors were calculated 
analytically thanks to the available closed-form solutions for certain (2D and 3D) elastostatic exterior 
problems - e.g. those for a penny-shaped crack, circular hole, and spherical inclusion in an infinite solid. 


To the authors’ knowledge, however, analytical solution to (18) is unavailable. As a result, numerical 


evaluation of {VY^ and thus A is pursued within a BIE framework |34l[TT]. Consistent with the assumed 
dimensional platform, the computation is effected assuming i) penny-shaped fracture of unit radius A = l, 
ii) elastic medium with unit shear modulus and mass density (// = 1, p = 1), and iii) various combinations 
between the trial fracture parameters and the Poisson’s ratio n of the background elastic solid. 

In this setting, general known properties of polarization tensors are deployed along with optimization 
techniques to decipher the numerical results into a mathematical expression for A which will be used 
later to speculate the fracture boundary condition from the TS distribution. 

To solve (18) for given Ttnai G A BIE computational platform is developed on the basis of the 


regularized traetion boundary integral equation m where the featured weakly-singular integrals are 
evaluated via the mapping techniques in [34|. Without loss of generality, it is assumed that the origin 


8 









of ^ coincides with the center of a penny-shaped fracture suface, and that n' = 63 (see Fig. B.ll). 
A detailed account of the adopted BIE framework, including the regularization and parametrization 


specifics, is provided in Appendix B 


To validate the computational developments. Fig. |^a) compares the numerically-obtained nontrivial 
components of with their analytical counterparts along the line of symmetry in Ftnai assuming 

traction-free interfacial conditions {i^n = Ks = ^) and u = 0.35. Here is noted that, thanks to the 
problem symmetries, the variation of along the ^i-axis equals that of |V 2 F^ along the ^ 2 -axis. 

To illustrate the influence of iTtnai oa the result. Fig. |^b) shows the effect of shear specific stiffness 
Ks (assuming = 0 ) on the tangential fracture opening displacement a similar behavior is also 

observed concerning the effect of oa the normal opening component iVsP^. 


Structure of the polarization tensor. Before proceeding further, it is useful to observe from (19) that only 
the part of A with minor symmetries enters the computation of T thanks to the symmetry of af and a. 
Further, as shown in nun], properties of the effective polarization tensor (hereon denoted by can 
be extended to include the major symmetry. With reference to the local basis ( 01 , 61,63 = n'), on the 


other hand, one finds that: i) =0 (a, /3 = 1 , 2 ) due to a trivial forcing term in (18); ii) oc 63 

owing to the symmetry (about the ^3 = 0 plane) of the boundary value problem for a penny-shaped 
fracture in solved by (18), and hi) 63 = 63 = 0 due to the anti-symmetry of the germane 

boundary value problem about the f 3 = 0 plane - combined with the axial symmetry of Ftnai about the 
f 3 -axis. A substitution of these findings immediateiy verifies that and are independent 

of whereas does not depend on Kg. Note, however, that the above arguments are predicated 

upon the diagonai structure of iTtnai according to (§. As a resuit, one finds that the effective poiarization 


tensor, superseding A in (19), permits representation 


A®" = as{K.s,p)'^{e3<Siep + ep®e3) <Si {e3®ep + ep®e3) + an{Kn,i^)ie3®e3®e3®e3). (20) 

p=i 

To evaiuate the dependency of ag and on their arguments, A is evaluated numerically according 


to (19) for various triplets {y^ng^Kn)- On deploying the Matlab optimization toolbox, the coefficients 
as{Ks, u) and an{Kn: ^.re found to be rational functions of their arguments, identified as 


ag{ng,v) = 


4(1-z/2) 


3(2 — n^i^Kgh. + z/ + 1) 


(^n: ^) — 


8 (l-z^)( 2 z/ + l) 
T + 1) 


( 21 ) 


where the impiicit scaling parameter, A = l, is retained to facilitate the forthcoming application of ( 21 ) 


to penny-shaped fractures of non-unit radius. Assuming n = 0.35, the behavior of ag and an according 
to (21) is plotted in Fig. |^c) versus the germane specific stiffness, with the corresponding numerical 
values included as dots. For completeness, a comparison between and the BIE-evaluated values 
of ag and an is provided in Eig. [^for a range Poisson’s ratios, n G [0.05,0.45]. 


Erom (20), it is seen that the interfacial condition on Ftnai has no major effect on the structure 
of the polarization tensor. This may explain the observation from numerical experiments (Section 


that by using a traction-free trial crack {Kg = Kn = 0 ) in ( 20 ), the geometrical characteristics of a 
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(a) ^ (b) (c) 








Figure 2: Dependence of the crack opening displacement and coefficients of the effective polarization tensor on 

the specific stiffnesses, Ks and /Cn, of a trial fracture assuming u = 0.35: (a) analytical values (solid lines) vs. numerical 
values (dots) of the COD along (^i,0, 0) for a traction-free crack; (b) evolution of the shear COD iVip^ with increasing 
ks, and (c) proposed variation (solid lines) vs. numerical variation (dots) of an and an versus the relevant specific stiffness. 


hidden fracture - its location, normal vector and, in the case of high excitation frequencies, its shape 
- can be reconstructed regardless of the assumed interfacial condition. Nonetheless, a trial crack with 
interacting surfaces introduces two new parameters {hZs and i^n) to the reconstruction scheme, whereby 
further information on the hidden fracture’s interfacial condition may be extracted. In this vein, it can 
be shown that the first-order topological sensitivity ([T^ is, for given a monotonic function of 
accordingly, precise contact conditions on F (the true fracture) cannot be identified separately via e.g. a 
TS-based minimization procedure. In principle, such information could be retrieved by pursuing a higher- 
order TS scheme [e.g. [12] which is beyond the scope of this study. As shown in the sequel, however, the 
present (first-order) TS sensing framework is capable of qualitatively identifying the ratio between the 
specific shear and normal stiffnesses on F - an item that is of great interest in e.g. hydraulic fracturing 
applications [e.g H HZ]. 



Figure 3: Performance of the formulas ( |21| describing an{i^n^i^) and as{i^s^i^) in the (i/,/^s,/^n)-space: numerical values 
(dots) vs. proposed expressions (solid lines). 


4. Qualitative identification of the fracture’s interfacial condition 


In this section, an ability of the TS indicator function (19) to qualitatively characterize the interfacial 
condition of nearly-planar fraetures is investigated in the “low” frequency regime, where the illuminating 
wavelength exceeds the length scales spanned by F. Such limitations are introduced to facilitate the 
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asymptotic analysis where the true (finite) fracture is approximated as being penny-shaped. The pro¬ 
posed analysis naturally extends to arbitrary-shaped, near-flat fractures by assuming a diagonal, to the 
leading order^ interfacial stiffness matrix and the COD profile proposed in m to calculate the relevant 
polarization tensor required in the asymptotic approximation of the scattered field due to T. In the 
approach, it is also assumed that the fracture location and normal vector to its surface are identified 
via an initial TS reconstruction performed using a traction-free trial crack {ks = n,n = 0) described 
in Section As shown there, however, the latter geometric reconstruction is not limited to a particular 
frequency range. 

On recalling the least-squares distance function (p in and the leading-order perturbation of J(0) 
in the TS may be rewritten as 


e^T(^”;n',K„,Ks) =-Re 


JSohs 


( 22 ) 


where (•)* signifies complex conjugation and —noting that u = u — 
is the scattered field due to T, see 0 . Given the fact that the hidden fracture is separated from the 
observation surface, the scattered field on 5'°^® can be expressed via a displacement boundary integral 
representation as 

■u(4) = y (ItiK®) 0n): dS'a,, ^ e (23) 

Small eraek asymptoties. Consider the testing configuration as in Fig. and let T with interfacial 
stiffness ^ be illuminated by a “low-frequency” plane wave in that kL ^ 1, where 2L is 

the characteristic size of T. With such premise and hypothesis that 2L < Ri ^ R 2 made earlier (recall 
that R 2 is the radius of 5'°^®), the hidden fracture can - in situations of predominantly flat geometry and 


an 0(1) aspect ratio - be approximated as a penny-shaped fracture of finite radius L. In this setting (23) 
can be expanded, utilizing the developments from Section as 


^(0 








4 e (24) 


where a (to be determined) is an indicator of the fracture location, and is its effective (low-frequency) 
polarization tensor given by 


AF = ^(^5(71(8)6/3-h 6/3(8) n) (g) (n(8)e/3-h 6/3(8) n) -h ^n(^G)n(8)n(8)n). ( 25 ) 

/ 3=1 

Here n is the normal on T; ( 61 , 62 , n) make an orthonormal basis, and 




4 ( 1 -z/ 2 ) 


3(2 — v^ikgL + z/ + 1) 


^n{kn 5 ^) — 


8 (l-z/)( 2 z/ + l) 
3{knL + 2z/ + 1) 


(26) 


where ks and kn are the shear and normal specific stiffness of the true fracture according to ([^, see also 
0 and ( [^ . On taking without loss of generality ( 61 , 62 , n) as the basis of the global coordinate system. 
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(24) can be rewritten in component form as 


Uj($) = ^ G (27) 

where /3 = 1 , 2 and the summation is assumed over repeated indexes as before. 

Assuming that the location and (average) normal on F are identified beforehand, one may set n' = n 
and ^° = 2 ; in ( 22 ) and expand the scattered field due to vanishing trial fracture F^ = + ^F^Hai 

analogous fashion as 

^i(C) = [4as{Ks,ty)T,lp{^,z,Lo)ali3{z)+aniKn,ty)^i3{^,z,uj)a(3{z)], ^ G (28) 


On substituting (27) and (28) into ( 22 ), the leading-order TS contribution in the low-frequency regime 
is obtained as 


T{z,Kn,Ks) ~ -L^ [an<;nQl{z) + (q!s + Q!n q) Q 2 (^) + asQQ 3 (^)], 


(29) 


where 


Qiiz) = l4wp/ [^il^i3n,z,u;)dS^, 

J Sohs 

Q2{z) = 

{ [4: 4 ] (•") 4 ] , 


(30) 


Qsiz) = 16 Re 


where (•)* indicates complex conjugation, a, /3 = 1 , 2 , and j = 1, 2, 3 as before. With reference to |Appendix A[ 
integration of the anti-linear forms, featuring components of the fundamental stress tensor, over 5'°^® can 
be performed analytically by approximating the distance r = |^ — zj, ^ G 5'°^® in (A.l) and (A. 2 ) to the 
leading order as r = ^ 2 - As a result, the behavior of (30) can be approximated as 


Qi — 


^33 




47rR2 


iX2 + |Re(X3) + |P(R2)|'], Q2 = 0(4) q 


^3 + I |F(R2)|"1, (31) 


owing to the initial premise of “far field” sensing (namely R 2 ^ 1 within the adopted dimensional 
platform), where F(-), P('), X 2 and X 3 are given by ( |A. 2 [ ) and (A. 6 ) in [Appendix A| along with the 
details of the calculation procedure. 


A remarkable outcome of the analysis is that the coefficient Q 2 describing the mixed term in (31) 
vanishes to the leading order, whereby the TS structure decouples and may be perceived as a superposition 
of the normal and shear contributions. Specifically, (29) becomes 


T{z;n,Kn,Ks) ~ {an{Kn,iy)<;nQi{z) + as{Ks,l')<;sQ3{z)), 


(32) 


where (^n and are given by (26). On account of (32) and the limiting behavior of an (resp. o^g) as a 
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function of the trial interface parameter Kn (reps. Kg) in (21), the ratio between the shear and normal 
specific stiffness (kg/kn) of a hidden fracture can be qualitatively identified via the following procedure. 


Interface characterization scheme. With reference to the last paragraph of Sectionit is hereon assumed 
that the fracture location is identified beforehand from the low-frequency scattered field data as 


= arg min^oT(C;n'(^°),0,0) 


(33) 


where, for each sampling point, is the optimal unit normal (minimizing T at that point) as described 

in [To]. In this setting the effective normal to a hidden fracture, n, is obtained by averaging n'(^°) over a 
suitable neighborhood of 2 ;. Here it is for completeness noted that the spatial variation of n'(^°) provides 
a clue whether the hidden fracture is nearly-planar and thus amenable to the proposed treatment, see 
the low-frequency results in Figs. and as an example. 

In the vicinity of 2 :, the TS characterization of the fracture’s interfacial condition is performed us¬ 
ing two (vanishing) trial fractures: one allowing for the COD in the normal direction only by assum¬ 
ing {Kg^Kn) = (oc,0), and the other restricting the COD to tangential directions via (Kg^Kn) = (0,oo). 


The resulting TS fields are then normalized by the relevant components of crf{z) according to (30). 
Thanks to the limits lim^^^^oo = 0 and lim^^^^oo = 0, this leads to the respective indicator func¬ 
tionals 


Ti(r) = , f c. -^a„(0,Of|X2 + |Re(X 3 ) + |P(i?2)P 


0-33 (^) 


T2(r) = 


4 r3^ r 

T(r;n,oo,0) ~-+ | |F(i? 2 )| 


/* / w X ■ VS 5 — d 2 ^sV^5 15 ^^21 3 I 


(34) 


were is in a neighborhood of 2 :. In practical terms, the latter is identified as a ball (centered at z) 
whose radius is a fraction of the germane (compressional or shear) wavelength. On the basis of ( [M| ), one 
can identify three distinct interface scenarios: 


1. The situation where kg and kn are of the same order of magnitude (e.g. traction-free crack), in which 
case [Ti/T 2 ](^°) = 0(1) regardless of u. This is verified in Fig. where the ratio T 1 /T 2 is plotted 
against R 2 - scaled by the shear wavelength = 27rCgfuj - for various Poisson’s ratios and two 
“extreme” sets of interfacial stiffnesses, namely k^ = 0.1 and k^ = 10 (^ = s,n). In the context 
of the proposed characterization scheme. Fig. [^a) plots the spatial distribution of Ti and T 2 in a 
neighborhood of a traction-free fracture. As can be seen from the display, the fracture is visible 
from both panels, which suggests that kg and kn are comparable in magnitude. 


2. The limiting case kg^kn that, in the context of energy applications, corresponds to a hydraulically- 
isolated fracture [HIITIIT]. Under such circumstances one has Ti(^°) <C T 2 (^°), which can be 
verified by letting ^ 0 in (26). This behavior is shown in Fig. [^b), where a hidden fracture 
with {kg^kn) = (2,100) is visible in the distribution of T 2 , but not in that of Ti. Note that the 
image of a fracture is notably smeared due to the use of low-frequency excitation, as postulated by 
the interface characterization scheme. 
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3. The case when the fracture surface is under small normal pressure and the effect of surface roughness 
is significant, namely kg^kn- Here Ti(^°) ^ T 2 (^°) due to the fact that lim/g^^oo = 0 thanks 


to (26). This is illustrated in Fig. [^c), where the fracture with (kg^kn) = (100,2) appears in the 
distribution of Ti only, suggesting that its interfacial stiffness according to (§ is dominated by the 
shear component kg. 


Here it is worth noting that the above TS scheme for qualitative identification of the ratio kg/kn is 
non-iterative, and shines light on an important contact parameter at virtually no computational surcharge 
- beyond the effort needed to image the fracture. In particular, since the (low-frequency) free and adjoint 
fields are precomputed toward initial estimation of the fracture location 2 ; and effective normal vector n, 
they can re-used to compute Ti and T 2 via (19), (20), (21) and (34), wherein the only variable is the 
effective polarization tensor - describing trial fractures with different interfacial condition. 



Figure 4: Ratio \Ji/~^ 2 ]{z) vs. the radius of measured in shear wavelengths, assuming the fracture specific stiffnesses 

as (a) ks = kn = 0.1, and (b) kg = kn = 10. 


Illumination by multiple incident waves. In many situations the imaging ability of a TS indicator func¬ 
tional can be improved by deploying multiple illuminating wavefields, which in the context of this study 
translates into multiple directions d of plane-wave incidence. In this case the “fortified” TS functional 
can be written as 


t = / T|d-w(d)d5d, 

J Qd 

which superimposes (in a weighted fashion) the TS distributions for incident plane waves spanning a 
given subset, of Hie unit sphere. In the context of (29) and (30), the only d-dependent items are the 
components of the free-field stress tensor ay. As a result, the criteria deduced from (34) remain valid 
under the premise of multiple incident-wave illumination provided that the free-held terms 


r /* / 1 
W33 ^33]^ 



a: 


/* 

3a 


^3/3 


in (30) are replaced respectively by 



^33] d 



^ 3/3 


]d 



^3/3. 
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5. Numerical results 


A set of numerical experiments is devised to illustrate the performance of the TS for elastic-wave imaging 
of fractures and qualitative characterization of their interfacial condition. To this end the BIE com¬ 
putational platform, described in Appendix B, is used to generate the synthetic data for the 

inverse problem. While the main focus of the study is on the low-frequency sensing as mandated by the 
characterization scheme, the results also include the TS reconstruction examples at intermediate-to-high 
frequencies, motivating future research in this area. The sensing setup, reflecting the assumptions made 
in Section 2, is shown in Fig. where the “true” fracture T is either i) a penny-shaped crack with di¬ 
ameter L = 0.7 and normal n = (0, —1/\/2, l/\/2) - see Fig. |^a), or ii) a cylindrical crack of length 
L = 0.7 and radius R = 0.35 shown in Fig. |^b). The shear modulus, mass density, and Poisson’s ratio 
of the background medium are taken as /i = 1, p = 1 and z/ = 0.35, whereby the shear and compressional 
wave speeds read Cg = 1 and Cp = 2.08, respectively. The elastodynamic field produced by the action of 
illuminating (P- or S-) plane wave, propagating in direction d, on F is measured over 5'°^® - taken as the 
surface of a cube with side 3.5 centered at the origin. This was done to investigate the robustness of the 
interface characterization scheme developed in Section]^ with respect to the simplifying assumptions used 
to derive (29)-(31), namely the premise that 5'°^® is a sphere. On the adopted observation surface, the 
density of sensing points is chosen to ensure at least four sensors per shear wavelength. In what follows, 
the TS is computed inside a sampling cube of side 2 centered at the origin; its spatial distribution is 
plotted either in three dimensions, or in the mid-section Hi (resp. 112 ) of tho penny-shaped (resp. cylin¬ 
drical) fracture shown in Fig.[^ In the spirit of an effort to test the robustness of the adopted simplifying 
assumptions, it is noted that the ratio between the radii of spheres circumscribing the sensing area and 
the sampling region is R 2 /R 1 = 1-75 ^ 1. 



Figure 5: Model problem: (a) sensing configuration with an embedded penny shaped fracture, and (b) with non-planar 
scatterer i.e. cylindrical fracture. 


Low-frequency TS reconstruction. Consider first the case of an “isolated fluid-filled” planar fracture 
shown in Fig. [^a), whose specific stiffnesses are given by (kg^kn) = (2,100). To geometrically iden¬ 
tify the fracture, the region of interest is illuminated by twelve P- and S- incident waves propagat¬ 
ing in directions d G {(± 1 , 0 , 0 ), ( 0 , ± 1 , 0 ), ( 0 , 0 , ± 1 )}, while assuming = (O 5 O) foi* fhe (van¬ 

ishing) trial fracture. The illuminating frequency is taken as cj = 4, whereby the ratio between the 
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probing shear wavelength = ^'kCqIu ^1.6 and fracture diameter, L, is approximately 2.2. For each 
incident wave and each sampling point, the optimal normal vector is estimated by evaluating 

T(r;n' = [cos(^i) cos(0j), sin(^^) cos(0j), sin(0j)], 0, 0) for trial pairs {0i^(j)j) G [0, 27r] x [0,7r] and choos¬ 
ing the pair that minimizes T(^°; •). Thus obtained TS distributions are superimposed as 


12 

t(r;n',0,0) = |min^ot|-i ^T(r;n', 0 , 0 )|(b„,d„), 

n=l 

resulting in a composite indicator function whose spatial distribution shown in Fig.j^a). In this setting, 
the fracture is geometrically identified via a region where T attains its most pronounced negative values, 
see Fig. ©b). To provide a more complete insight into the performance of the approach. Fig. |^c) plots 
the distribution of point-optimal unit normal over the reconstructed region. For the purposes 

of interface characterization, the fracture location 2 : is identified according to ( [^ , while its effective 
normal n' is computed by averaging over the reconstructed volume in Fig. [^b). 



Figure 6: Low-frequency reconstruction of an “isolated fluid-filled” planar fracture with (kg^kn) = (2,100): (a) composite 
TS distribution T(^°; n', 0, 0) in three dimensions, and in the fracture mid-section Hi; (b) region containing the most 
pronounced negative TS values: —1 ^T^ —0.6, and (c) point-optimal normal vector plotted over the true fracture 

surface F. 


For completeness, the above-described geometrical identification procedure is next applied to the 
cylindrical fracture in Fig. [^b), assuming the “true” specific stiffnesses as (kg^kn) = (4,4) and setting 
the excitation frequency to cj = 5. In this case, the wavelength-to-fracture-size ratio can be computed 
as Xg/L ^ 1.8. The resulting distributions of T and point-optimal normal n' are shown in Fig. from 
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which one can observe that i) the method performs similarly for both planar an non-planar fractures, 
and ii) the reconstruction procedure is apparently not sensitive to the nature of the interfacial condition 
in terms of kg and kn- 



Figure 7: Low-frequency reconstruction of a cylindrical fracture with (kg^kn) = (4,4): (a) composite TS distribution 
T(^°;n',0, 0) in three dimensions, and in the fracture mid-section 112 ; (b) region containing the most pronounced negative 
TS values: —1 —0.6, and (c) point-optimal normal vector plotted over the true fracture surface F. 


Interface characterization. With reference to Fig. [^a), the characterization of a penny-shaped fracture 
with three distinct interface scenarios is considered, namely: i) the traction-free crack i.e. = 0, 

ii) isolated fluid-filled fracture with (/cg, kn) = (2,100), and iii) fracture with rough surfaces under insignif¬ 
icant normal stress [42], simulated by setting (/cg, /c^) = (100, 2). In all cases, the interface characterization 
is carried out at u; = 5 i.e. Xg/L 1.8 using twelve incident (P- and S-) waves as described earlier. Next, 
to identify the contact condition, the “test” indicator functions 


ti(r) = in 




E 12 

n=l 


^33 




t2(r) = in 




^12 r / 

\(7ia(T. 


Fill [c 




^ ] W| 


(hn,dn) 


are computed on the basis of (34), where Tm = min{min^oTi, min^oT 2 }. The resulting distributions of Ti 
and T 2 for all three scenarios are plotted in Fig. using a common color scale, over the fracture’s mid¬ 
section Hi. From the display, it is clear that Ti(|°)/f 2 (|°) = 0(1), Ti(C°) ^ 2 ( 1 °), and Ti(|°) > t2(C) 
respectively for the “true” interfacial conditions according to i), ii) and iii). These results indeed support 
the claim of the preliminary characterization scheme that, at long illumination wavelengths, the interfacial 
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condition of nearly-planar fractures can be qualitatively assessed at virtually no computational cost - 
beyond what is needed to identify the fracture geometrically. 


TS reconstruction at higher frequencies. Clearly, the geometrical information in Figs, [^b) and [^b), 
obtained with Xg/L ^ 2, does not carry sufficient detail to accurately reconstruct the fracture surface 
in 3D from the scattered elastic waves. To help mitigate the drawback, the forward scattering problem 
is recomputed at a higher frequency, namely oo = 20 for which Xg/L 0.45. In doing so, the number of 
incident elastodynamic fields is increased so that twenty plane waves of each (P- and S-) type, propagating 
in directions n=l,2 ,...,20}- evenly distributed over the unit sphere - participate in the TS 

evaluation. The resulting “high-frequency” behavior of the composite indicator function 

40 

t(r;n',0,0) = |min^ot|-i ^T(r;n',0,0)|(5„,d„) 

n=l 


is shown in Figs.|^and 10 for the cases of penny-shaped fracture and cylindrical fracture, respectively (see 
Fig.§. The featured results are consistent with the recent findings in acoustics and elastodynamics [221 
SSIIIT] which demonstrate that, at higher illumination frequencies, pronounced negative values of TS 
tend to localize in a narrow region tracing the boundary of a scatterer. However, it the present case T 
also exhibits a notable sensitivity to the fracture’s interfacial condition; in particular, for the traction-free 
crack in Fig.j^a), extreme negative values of the TS are localized in the vicinity of the crack tip^ whereas 
in the case of fractures with interfacial stiffness - Fig.j^b) and Fig. 10 - the TS exposes the entire fracture 
surface. These initial results suggest that at shorter wavelengths, the TS experiences a different type of 
sensitivity to the fracture interfacial condition, which may lead to a more detailed identification of the 



Figure 8: Spatial distribution of Ti (top panels) versus T 2 (bottom panels), in the mid-section of a penny-shaped fracture, 
whose interface is (a) traction-free (kg = kn = 0), (b) isolated fluid-filled (kg = 2, kn= 100), and (c) of significant surface 
roughness and under low normal stress (kg = 100, kn = 2). 


18 
















Figure 9: “High”-frequency (Ag/L ~ 0.45) TS reconstruction of a penny-shaped fracture whose interfacial condition is (a) 
traction-free (kg = kn =0), and (b) isolated fluid-filled (kg = 2, kn = 100). 


fracture’s specific stiffnesses. Given a sensory data set that includes the scattered field measurements 
at both long and short wavelengths, one may consider a staggered approach where i) high-frequency 
data are deployed to more precisely evaluate the fracture geometry, including its location 2 : and effective 
normal n; ii) low-frequency observations are used as in Section]^ to qualitatively identify the interfacial 
condition, and iii) additional information is obtained on the fracture’s interface thanks to the dependence 
of the high-frequency TS thereon. Such developments, however, require high-frequency asymptotics of 
the scattered field due to a fracture with specific stiffness - a topic that is beyond the scope of this study. 



of T in the fracture’s mid-section, and (b) corresponding T map thresholded at 70%. 
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6. Summary 


This work investigates the utility of the topological sensitivity (TS) approach as a non-iterative tool for 
the waveform tomography of fissures with specific stiffness, e.g. hydraulic fractures. On postulating the 
nucleation of an infinitesimal penny-shaped fracture with constant (normal and shear) interfacial stiff¬ 
nesses at a sampling point, the TS formula and affiliated elastic polarization tensor are calculated and 
expressed in closed form. In this setting, it is shown via numerical simulations that the TS carries the 
capacity of exposing the fracture location and its unit normal from the long-wavelength scattered field 
measurements, regardless of the assumption on the (trial) interfacial parameters of a vanishing fracture. 
Given thus obtained geometric information and elastodynamic (free- and adjoint-field) simulations re¬ 
quired for its computation, it is further shown that the interfacial condition of nearly-planar fractures can 
be qualitatively identified at virtually no added computational cost, using two auxiliary TS maps evalu¬ 
ated for certain (extreme) combinations of the trial contact parameters. In particular, the analysis shows 
that such scheme allows for the ratio between the shear and normal specific stiffness - representative of a 
hidden fracture - to be exposed as either i) near-zero, ii) on the order of unity, or iii) exceeding unity by 
a large amount. Such information can be used both directly, e.g. to discriminate between the old, newly 
created and proppant-filled fractures, and as an “initial guess” on the fracture’s interfacial condition for 
a more comprehensive, nonlinear optimization approach to waveform tomography. Through preliminary 
simulations at “short” incident wavelengths - subpar to the characteristic fracture size - which demon¬ 
strate both enhanced imaging resolution of the TS indicator function and its heightened sensitivity to 
fracture’s interfacial condition, this study further provides an impetus for studying the high-frequency 
elastodynamic scattering by, and TS sensing of, fractures with specific stiffness. 
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Appendix A. Asymptotic behavior of the integrals in (30) 


The aim of this section is to expose the leading-order behavior of the integrals in (30) for R^jRx ^ 1, 
where the integrands are expressed as antilinear forms in terms of the components of the fundamental 
stress tensor. 


Elastodynamic fundamental solution. Assuming time-harmonic excitation by the unit point force applied 
at cc G in direction k, the governing equation of motion for an infinite solid with shear modulus /i, 
mass density p and Poisson’s ratio v can be written as 
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where U^{^,x,uj) and x,cj) denote respectively the fundamental displacement vector and stress 

tensor, given in the component form by 


Ui{i,x,Lo) = ■^^{Bi{r)5ik +B 2 {r)r^ir^k), 

= ^^(2B3(r)r,ir,jr,fc + {5ikrj + 5jkr,i)B4^{r) + B5{r)8ijr^k)- 


(A.l) 


Here 
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iXp)e ‘ 


-iXp 


(A. 2 ) 


dr 

Wi' 


where Sij denotes the Kronecker delta, Cg = y//i/p, and Cp = Cg y/{2 — 2u)/{1 — 2u). 

Integration. For x G and ^ G 5'°^® (see Fig. [^, the argument r = \^ — z\ of the relevant coefficients 
in (A.l) and (A. 2 ) describing can be approximated (assuming R^jRx ^ 1) to the leading 


order as r i? 2 ; as a result, the prevailing behavior of integrals in ( |3Q[ ) can be written as 

c p^iTT 

Qi-. / ~ Xi / / (X2rl + 4Re(X3)r| + |B5(i?2)P)sin(0)d0d,/., (A.3) 

Jo ^0 

r r^iTT 

Q2-. / [Ei3EL](^,;s,:^)d5^ ~ Xi / / r3r7X2r23 + 2X3)sin(,/.)d0d,/. = 0, (A.4) 

J Jo Jo 


5 'obs 


n ZTT 

(X 2 r >2 + (r| +r2)|B4(i?2)|2)sin(,/.)d0d,/., 


Q3: - (A.5) 

/ [SfiSy(4,;2,a;)d54 ci Xi / / (X 2 ri r 2 r| + 154(^2) |" rir 2 ) sin((/.)d0d(/. = 0, 

Js'obs Jo JO 

where /3 = 1 , 2 and 


Xi = 


r, X2 = 4Re{[|R31'+ 2^3*54] (i?2)}, X3 = [IR4P + (^3 + R4)B;](R2). (A.6) 


(47ri?2)2’ 

Considering the unit vector Vr = ( 71 , 725 ^, 3 ) = (sin(0) cos( 6 >), sin(0) sin( 6 >), cos(0)) used to define 5'°^® 
in spherical coordinates, the integrals of (|A.3[)-([A^ are analytically evaluated which results in (32). 
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Appendix B. BIE computational platform 


This section summarizes the numerical scheme adopted to solve the elastodynamic traction BIE for 
a fractured three-dimensional solid. The approach borrows substantially from the ideas established 
in considering the regularization of the featured surface integrals. For brevity, the technique is 

described with reference to the elastostatic canonical problem Q. With slight modifications, however, 
this method is utilized in a more general setting of Section to compute the scattering of elastic waves 
by an arbitrarily-shaped fracture. Accordingly, the auxiliary formulae are expressed in their most gen¬ 
eral form (as applicable). Unless stated otherwise, the Einstein summation convention is assumed over 
repeated coordinate indexes. 


Regularization. To avoid evaluating the Cauchy principal value in (18), the featured integral equation 
can be conveniently rewritten as 


-[ei0e,-+e,-®e4,n^ - (0 = + 


[ {D,,lV^r{x) - D,4K„F^(0) x) dSs., 


(B.l) 




where the dummy indexes are summed over 1,3, and the singularity is transferred to the auxiliary integrals 

comprising the third-order tensor 
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Here it is useful to recall that r = —x|, i/ = di/dxi, and = Kn{es0es)-\-Ks{ei0ei)-\-Ks{e2(S)e2) 

where (ei, 62 , 63 ) are the unit vectors along (fi, C 2 , fs)- Considering Ii first, one finds [TT] via integration 
by parts, taking advantage of the Stokes identity, and noting that = 0 on STtnai, that the first 


in (B.2) can be reduced as 


If = 


[ 7 d5s - / 

■/rtrial ^ ^ JdTt, 


ds, 
. r 


ler. 


(B.4) 


where Ttnai is interpreted as an open set (excluding fracture’s edge QTtHai); = 1,2,3; v denotes the 
outward normal on STtnai lying within the tangent plane to Ttnai, and Wk{f) = /,/c — '^kif^p'^'p) i^ the 
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tangential derivative operator. In terms of Jpqi, on the other hand, one similarly obtains 


Spq I \n'fn'^i^s dSs: + I - n'f{2npn'g - T^pT^q)wsn'^ dSg;- 

"'rtPial ^ .frtPial ^ 

2/ -^n'gn'in'gr^s dSx+ [ -(v^n^n'^ + v^n^n^ - r^^r,,v<;) ds + 

•/rtHal •fartHal^ 

I -'>T''e^,q^,s{vsnp - Vpn'J ds, I e rt^ai- 

Jdruui r 

(B.5) 

For further details on (B.1)-(B.5), the reader is referred to Chapter 13 in m- Here one should mention 


that, for the canonical problem in Section 3.2 where n' = 63 , formulas (B.4) and (B.5) can be remarkably 


simplified (see Appendix 5.A in m)- In this presentation, however, both formulae are kept in their 
general format for they also pertain to the scattering by arbitrarily-shaped fractures. 


Pammetrization. With reference to Fig. |B.ll| the fracture boundary FtHai is discretized using a conformal 
mesh permitting surface parametrization {y ^ x) as 


x{yi,y 2 ) = i^m{yi,y 2 ) m = l,-",Nn, 


(B. 6 ) 
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Figure B.ll: Geometric (dotes) and interpolation (crosses) nodes on Ftriai in physical and parametric spaces. 


Here = 8 designates the number of nodes per element, and denotes the global coordinates of 
the element’s m^^ node - whose shape function is that of the standard eight-node quadrilateral 

element. In this setting, one finds the natural basis ai ,2 of the tangent plane and the surface differential 
as 

= ^(y)d5y, 5^(y) = [aixaal, (B.7) 


, , -r, 

a/ 3 (y) = —X 


dy/3 

where (3 = 1^2 and the dummy index m is summed over 1 , 8 . At this stage, one should note that all inte¬ 
grands in (B.4) and (B.5) - comprising in (B.l) - are known so that the boundary parametrization. 
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given by (B. 6 ) and (B.7), is the only necessary step prior to numerical integration. 


In light of the smoothness requirement by the traction BIE (B.l) and the adverse presence of the 


tangential derivative operator Dqs{-)^ the COD |V]*-^(x) is discretized via non-conforming interpolation 
(see Section 3.2 of m for details). In particular, the interpolation i.e. collocation points are situated 


inside the boundary elements (see Fig. B.ll), and their position with respect to the geometrical nodes - 
in each element - is quantified via parameter 77 in the (^ 1 ,^ 2 ) space. In this setting, the COD over the 
parent element can be approximated as 


[F,F^‘(X) = My) mu. m = 1,...,8, 


(B.8) 


where 


WifA = = -f^{r] + yTyi)iji + y^y2){yfyi + y2y2-'n). (b.9) 


B.ll 


and (^i^,^ 2 ^) is the position of the m^^ collocation point in the parent element as shown in Fig. 

Here it is important to mention that the surface elements adjacent to ^FtHai ^.re of quarter-node type 
(see e.g. Chapter 13 of m), designed to reproduce the square-root behavior of the COD in the vicinity 
of QFtriai- Note that for constant distribution of interfacial stiffness, i.e. for constant stiffness matrix Xtriai? 
the asymptotic behavior of the COD near ShtHai remains the same as that in the case of traction-free 
crack (see [46] for proof). 


Given (B.6)-(B.9), it can be shown m that the tangential derivative operator DqslVM{x) in ( jB.lj ) 
permits the parameterization 

D.slVeFix) = A” (y) A“ (y) = ^ [{a2-ep)cPr.,i - (B.IO) 

where the dummy indexes p and m are summed over 1 ,3, and 1 , 8 , respectively; the basis unit vectors Cp 
are shown in Fig. B.ll, and Cpqs denotes the Levi-Civita symbol. 


On substituting (B.7) and (B.IO) into (B.l), one arrives at the algebraic system for the values of 
at the collocation nodes | as 

i[e,0e,.+e,.0e4,n^ - A4qy“*)54(5y,r,,,)[K,l4. + 

F F / (^qs(y) -'^ee-A™®(y'"*))E4(sy,S(y))^e(y)d4}l^mlme- (B.ll) 

e=l m=l ^ ' 

^(1 -4e.)E Y'iCkepf E?;,,(sy,*(y))5^e(y)d5jA-f(y-*)}[y„l4. , 

e=l m=l 

where e* is the element number; m* denotes the local node number; no sum over e* and m* is implied; 
Nn = 8 ; Ne is the number of elements, and 

ML = iVifHxf), A^:{y) = |g^[(a^ep)<).„,i-(4-ep)0„,,2](y). (B.12) 


Here it is worth noting that all integrals in (B.4), (B.5) and (B.ll) are numerically integrable. A specific 
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mapping technique (in the case of weak singularity) along with the standard Gaussian quadrature method 
is employed to evaluate the aforementioned integrals (see Section 3.9 of [24] for details). 
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